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Abstract 

This paper formed part of a preliminary research report for a risk consultancy and 
academic research. Stochastic Programming models provide a powerful paradigm for 
decision making under uncertainty. In these models the uncertainties are represented 
by a discrete scenario tree and the quality of the solutions obtained is governed by the 
quahty of the scenarios generated. We propose a new technique to generate scenarios 
based on Gaussian Mixture Hidden Markov Modelling. We show that our approach 
explicitly captures important time varying dynamics of stochastic processes (such as 
autoregression and jumps) as well as non-Gaussian distribution characteristics (such 
as skewness and kurtosis). Our scenario generation method enables richer robustness 
and scenario analysis through exploiting the tractable properties of Markov models 
and Gaussian mixture distributions. We demonstrate the benefits of our scenario gen- 
eration method by conducting numerical experiments on FTSE-100 data. 

Key words: Markov Processes, Risk Analysis, Stochastic Programming, Scenarios. 



1. Introduction 

Businesses are constantly confronted with the position of making decisions under 
uncert ainty e.g. portfolio investment, power production and even agricultural planning 
BLQ?! . Therefore the ability to be able to quantify and optimise such decisions is of 
paramount value and the modelling paradigm of Stochastic Programming offers such 
a possibility. Stochastic Programming (SP) has been proven to provide practical and 
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optimal solutions e.g. M ulvey's ce lebrated Towers Perrin-Tillinghast model saved US 
West $450-$ 1000 million 



MGMOO |. 



Stochastic Programming concerns decision making under uncertainty. A multistage 
stochastic programming problem is defined as: 



minEp{/(x, a;)} = min / f{x,uj)dP{uj), 



(1) 



where 



• x={xt} is the set of decisions at each stage t=l,2....T where T is the total number 
of stages; 

• S is the set of feasible decisions; 

• a; is the stochastic process; 

• f2 is the set of all possible events; 

• P is the probability measure for stochastic process u; 

• Ep is the expectation taken with respect to probability measure P; 

• /(x,a;) is the cost function. 

In stochastic programming a method is required to discretise the random process u 
into a set of discrete outcomes, known as scenarios. The method of obtaining scenarios 
is known as scenario generation. 

The root node of the scenario tree represents 'today' and is known with certainty. 
Nodes or scenarios at later stages progressing from the root node represent possible 
realisations of the rando m proc ess, where each scenario has a probability and a value. 



As Birge and Louveaux BL97l | point out, scenarios are particularly useful when the 
optimal solution varies considerably with changes in value of the stochastic variables. 
If this were not the case, the solution obtained for one stochastic variable realisation 
ought to be a good solution for most realisations. 

The general formulation of equation (1) under scenario generation is: 



i=Ni i=N2 



minEp{/(x,w)} = min V p\f{x\,uji) + V p^f{xluj'^) + 

1=1 

=Nt 



1=1 1=1 

i=NT 



T)-i 



1=1 

where 
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• i=l,2....Nt is the scenario and Nt is the total number of scenarios at stage t] 

• pI is probabihty of scenario i at stage t; 

• f{xl,ul) is the cost function of scenario i at stage t. 

A key issue in scenario generation that tends not to be addressed is the abihty to 
capture time varying dynamics of the random process uj (such as jumps and reversions) 
and to a lesser extent, capturing non-Gaussian characteristics. Our scenario generation 
method specifically addresses the aforementioned problems through the properties of 
Gaussian Mixtures (GM) and Hidden Markov Models (HMM). 

In GM HMM we model a time series as a HMM where each hidden Markov state 
is represented by a univariate mixture of Gaussians. To produce a scenario tree we let 
each Markov state correspond to a stage in the scenario tree. Then we obtain scenarios 
at each scenario tree stage by sampling from each state's associated univariate mixture 
of Gaussians. Hence the key part in our scenario generation procedure is determining 
the GM HMM and then selecting a state sequence for our scenario tree. Once the 
states have been selected, we simply sample from each state's GM to obtain scenarios 
at each stage. 



Guided Tour 

The outline of this paper is as follows: firstly we review current scenario generation 
methods and then introduce the GM HMM model. We then explain the advantages 
of the GM HMM for scenario generation, showing how the particular features of GM 
and HMM capture the time varying dynamics of data. We also show how GM HMM 
provides SP modelling benefits e.g. robustness analysis. Finally, we illustrate the 
GM HMM scenario generation and modelling benefits by applying it to the FTSE-100 
index, giving computational results. 



2. Current Scenario Generation Methods 



Currently there exist four main categories for scenario generation DDBMVl |: sam- 
pling, simula tion, stati stical property matching and hybrid methods e.g. si mulation- 



AGKL03| . For a detailed list of scenario generation methods see DDBMVl | , MitOGl ] 



optimisation 

In simulation, a stoch astic pr ocess's simulation results in a "path by path" produc- 
tion of the scenario tree KW03l |. To reduce the number o f paths sometimes paths are 
bounded together by some grouping method, for instance Igrso^ . Popular stochastic 
processes include geometric Brownian motion and its variants. 
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In scenario generation by sampling we sample values from a distribution and the 
sampled v alue represents the scenario's value. Various methods of sampling exist 



AGKL03l |. giving rise to different sampling scenario ge neratio n methods. Popular 



methods include: random sampling, importance sampling DT99[ |. bootstrap sampling, 
internal sampling KWOSj , conditional sampling and stratified sampling. 

In scenario generation by statistical property matching we do not require knowledge 
of the random variable's (w) pdf (probability density function). Instead we require the 
statistical properties of the generated scenario tree to match the statistical pr operties 
of some target (as closely as possible) e.g. percentiles. Hoyland and Wallace in [HWOll ] 
describe a scenario generation method using moment matching. 

3. Gaussian Mixture Hidden Markov Model Scenario Generation 

3.1. Markov Models 

Markov Models (MM) and Hidden Markov Models (HMM) are methods of math- 
ematically modelling time varying dynamics of some statistical process. They only 
require a weak set of assumptions yet provide powerful results. MM and HMM model 
a stochastic process (or any system) as a set of states with each state possessing a set 
of signals or observations. Changes in the stochastic process or system are modelled by 
moveme nts betw een states. The models have been use d in di verse applications such as 



mance 



MGPG06|. 



HarOl| . queuing theory SF06| . Engineering TGOl| and biological modelling 



A Markov model (or more specifically a Markov chain) is a stochastic process Xt 
with a countable set of states and possesses the Markov property. That is, given Xt is 
in state it (denoted Xt = it), Xt+i is not influenced by the set of values Xu for u < t. 
Formally, the Markov property is 



P{X^ 



j\Xo = io,Xi 



P{X, 



t+1 



j\Xt 



where it denotes a state at time t. As time passes the process may remain or change to 
another state (known as state transition). The state transition probability matrix (also 
known as the transition kernel) tells us these probabilities Pij = P{Xn+i = j \ Xn = i) 
where P^j is the probability of the process changing to state j given that we are now 
in state i. Note that we assume all probabilities are stationary in time. 

3.2. Gaussian Mixture Hidden Markov Models 

In Markov models we directly observe the state associated with each signal or 
observation; whereas in Hidden Markov models the states are directly unidentifiable 
or unobservable. However in HMM we assume that each signal or observation is a 
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probabilistic function of each hidden state, hence we can statistically infer the HMM 
states from observations. If observations were not statistical functions of state then we 
could not determine the hidden states. The observation probabilities are represented 
in the observation matrix B, showing each observation probability for each state. Note 
that as in plain Markov models we assume all probabilities are sta tionary. 
A hidden Markov model is fully defined when we know Rab89| : 

• A'' is the total number of hidden states Si, S2-.., Sj\f; 

• A represents the hidden state transition matrix of size NxN. Each element is 
ttij = P{qt+i = Sj\qt = Si) where qt is the hidden state at time t; 

• vTj (=P(go = Si)) for 1 < i < N denotes the the probability of b eing in state Si 
at time t=0; 

• C denotes the number of distinct observations. The set of observations are rep- 
resented as Vi,V2 ■ ■ ■ , Vc- For example if we model weather as a HMM with only 
2 observation signals hot and cold, then C = 2 and t'i=hot, f2=cold; 

• B denotes the observation matrix. Each entry is hj{c) = P{vc\qt = Sj), 1 < c < 
C. So bj{c) is the probability of observation Vc at time t, given that we are in 
state Sj at time t. 

A Gaussian mixture is a pdf composed of a weighted sum of univariate Gaussians, 
where all the weights sum to one. Let ^(<^) be a ra ndom variable ^{uj) G M whose pdf 
f{^) is represented by m Gaussian mixtures EH8l| : 

/(O = WlN{^, fii, (Ti) + W2N{^, fl2, (72) + .... + WmN{^, flm, (Tm) , (2) 

where fim, Wm denotes the mth Gaussian's mean, standard deviation and w eightin g 
respectively. Note that the mean and all other moments are weighted sums too EH8l| , 
for instance the overall distribution's mean fif is: 

m 

fJ-f = ^WkiXk- (3) 

k=l 

Now a Gaussian mixture HMM is simply a HMM where each state's observations is 
jy a G aussian mixture pdf. In other words fj{C} is our observation signal 



represented 
hj{c) where [LevOSl ] : 

m 

fi (0 = WjkN{i, fijk, CTjk) , e e M, (4) 

k=l 

where 
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• is the pdf of ^ in state j where I < j < N; 

• Wjk is the mixture weighting in state j for Gaussian mixture k where 1 < k < m, 
Wjk > Wj, k and YlT=i ^jfc=l Vj; 

• fijk is the mean for each mixture k, in state j; 

• ajk is the standard deviation for each Gaussian k in state j; 

• m is the total number of Gaussians in the mixture. 

3.3. GM HMM Scenario Generation Method 

Given we are applying GM HMM to financial data we may expect the stock market 
to belong to 3 economic phases (growth, recession and a transitional state) so we set 
N=3. Additionally, we choose m=2 given that most empirical financial return distri- 
butions can be described by a mixture of 2 Gaussians. We should also note by the 
Akaike information criterion that it is better to minimise the number of variables to 
prevent overfitting. 

Scenario Generation Algorithm 

Let f where f ~ t^[0, 1] denote a random number drawn from a unit interval. 

1. Randomly select an initial state: 
Draw a random number f. 
If: 

• f < TTi set j = l; 

• TTi < f < TTi + set j=2; 

• otherwise set j=3. 



2. Randomly select a Gaussian component k ( EH8l| . BDM7l| ): 
Draw a random number f. 

Set k=l if f < Wji, otherwise k=2. 

3. Sample Gaussian A^(yUjfc, cTjfc)- 

4. Randomly select the next state: 
Draw a random number f. 

If: 

• f < aji set j'=l; 

• flji < f < flji + aj2 set j'=2; 
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• otherwise set j'=3. 

Set current state to new state: j i-^ /. 
5. Increment I. 

If I equals the maximum number of scenarios required then stop, otherwise goto 
step 2. 

To calibrate the GM HMM we apply the Baum- Welch algorithm, which is an 
expectation-maximisation algorithm. The algorithm determines the GM HMM model 
parameters such that it maximises the statistical likelihood of a given observation data. 
The algorithm is guaranteed to converge to some local optimum f or calcu l ating the like- 



lihood of a given observation sequence. The reader is referred to Rab89l |. JR9l| . JR85[ | 
for more information. 

If we define A={A, B, tt,;} to denote the set of HMM parameters we wish to estimate 
then the pseudo code for the Baum- Welch algorithm is: 

1. Initialisation: 

Input initial values of A (otherwise randomly initialise) and calculate P{0\X) 
using the forward algorithm. 

2. Estimate new values of A. 
Iterate until convergence: 

(a) Using current A calculate jtihj) = P{<lt = Si,qt+i = Sj\0,X) for each t,i,j 
where 



at{i)aijbj{Ot+i)pt+i{j) 



Variables at(i), Pt(i) a re known as the forward and backward variables re- 



spectively (see |Rab89|): 



aS) = P{0i02 Ot,qt = S,\X), (6) 

PtU) = P{0t+i0t+2 OT\qt = S,,X). (7) 



(b) Calculate new A parameter estimates using 7t(i, j) (see |Rab89l | ) . 

(c) Calculate P{0\X) with new A values using the forward algorithm. 

3. Goto step 4 if two consecutive calculations of P{0\X) are equal (or converge 
within a specified range). Otherwise repeat iterations: goto step 2. 

4. Output A. 
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3.4- Scenario Generation Salient Features 

Explicit Modelling of Reversionary and Jump Dynamics 

Stochastic processes or time series tend to exhibit reversionary characteristics (e.g. 
autoregression) and jumps. These dynamics make continuous time modelhng challeng- 
ing and are typically ignored in scenario generation since most methods are unable 
to accommodate these dynamics. Using HMM we can explicitly capture reversionary 



characteristics and 



u mps; in fact HM M have been used widely in Finance to capture 



such aspects jRebOJ], HarOl[ |. HamSQj. Important reversion properties captured are: 



autoregression: this is captured by the transition probabil it ies of returning to the 
same current Markov state P{qt+i = Si\qt = Si) ( HamSol . HanJ |). 



conditional heteroscedasticity:since each Markov state has a different variance 
and each future state's transition depends on the current state, variance changing 
is conditional on the past. 



• cyclical behaviour [Ham89| : the cyclical nature of data is captured through the 
transition probability matrix A, so that the MM cycles through its entire set of 
states. 

Jumps are explicitly modelled by the transiti on prob ability of switching between 2 
states with significantly different distributions RebO^ . 



Modelling Diverse Distributional Characteristics 



■act GM 



TSMSSj . 



Gaussian mixtures can capture virtually any empirical distribution. In 
have been used to accurately capture outliers and "noisiness" in empirical data 
Various distribution characteristics are captured by simply adjusting the mixture com- 
ponent means, variances and weighting. 

Tractable Statistical Calculations 

Gaussian distributions tend to be apphed to data modelling due to their analytical 
tractability yet most (empirical) distributions tend to be non-Gaussian e.g. asymmetric 
and fat tailed. Additionally non-Gaussian models tend to possess intractable statistical 
calculations e.g. the Levy distribution. 

Using a mixture of Gaussians we can model most non-Gaussian distributional char- 
acteristics (e.g. skewness, kurtosis etc..) without losing tractability. For example, the 
cumulative distribution function for a GM mixture consisting of mean components 
is simply a weighted sum of each component's cumulative distribution value. Thus the 
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the overall cumulative distribution function for any univariate GM -F(^) gives BSOTj]: 

- tJ-k 



k=l ^ 



(8) 



where ^(z) is the standard normal cumulative distribution function for random variable 



z. 



Copula Approximation 

Copulas map marginal distributions to 1 multivariate distribution and so are par- 
ticularly useful in risk analysis. However copulas can become analytically intractable, 



causi ng calcu 
tion 



ation and computation difficulties e.g. sampling and scenario genera- 
HKW03l | . Furthermore, there exist significant implementation problems since 
we would need to l)estimate the covariance matrix, which tends to be unstable over 
time 2)specify the copula type (e.g. median, Gaussian) 3)determine a feasible method 
of scenario generation from the overall multivariate distribution. All these tasks are 
non-trivial. 



Mixture of Gaussians can be used as an alternative to copulas [BS07|, where all 
the component distributions are mapped into one univariate Gaussian mixture. Using 
univariate GM the problematic steps of 1-2 are removed, rather only the weighting for 
each component distribution need be estimated. 

If we map to a univariate (rather than a multivariate) distribution then the distri- 
bution model is more amenable to analytical calculations, computational methods and 
scenario generat ion. For instance GM can be applied to models of geometric Brownian 



motion 



BMS03l | whereas for multivariate Gaussians this becomes non-trivial. Univari- 



ate GM can be applied to a wider and more efficient range of scenario generation and 
computational methods. 

Robustness and Scenario Analysis 

In SP we m ay wish to evaluate a SP model's robustness by stress testing an d shock 
analyses |BL97l |. Dupacova suggests a contamination method for stress testing [Dup96l ] 
where randomness is added to the scenarios. For shock analysis we typically expose 
the SP model to an extreme event or a worst case scenario. 

Using GM we can enrich stress and shock analy sis, in fa ct GM already have been 
used as robustness methods in statistical modelling TSM85| but not in SP. For stress 
testing, we can add a GM of any choice to the overall distribution, representing a 
contamination as in Dupacova's method. Furthermore, we can calculate the impact on 
the overall distribution's moments and pdf for given weights, means and variances. 
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For shock analyses, rather than simply testing with respect to an extreme event 
scenario we can model a shock by a GM of chosen weighting, mean and variance 
and this is more realistic. For instance an overall GM distribution may represent a 
company's profit and we may wish to determine the impact of some extreme event 
e.g. hurricanes can severely disrupt natural gas supplies. We could more realistically 
model a hurricane event by a GM compared to testing with respect to one extreme 
event scenario. 

Adaptability to Multivariate Distributions 



The GM HMM can model multivariate Gaussian mixtures (see LevOSl ] for a detailed 
application in Engineering) and there exist significant statistical literature, algorithms 
and programs available. Additionally, multivariate GM HMM captur e time varying 
covariance matrices by specifying one for each hidden Markov state Rab89l ] . This 
is a significant advantage over other multivariate modelling methods particularly for 
financial applications, where correlations are critical to portfolio diversification. 



4. Numerical Experiment: Case Study Application to the FTSE-100 

We present results of applying the G M HMM to modelling the monthly returns of 
the FTSE-100 vt, where is defined by HulOOj ] : 

It 



rt = In 



(9) 



t-i 



The variable It is defined by the FTSE-100 index value in month t. We used monthly 
data from April 19 84-February 2007; the GM HMM was calibrated by the Baum- Welch 
algorithm Rab89l | using an implementation by Murphy MurOSj. The results are now 
presented. 
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Table 1: Transition Matrix A Results 



Current State 


Next State 


State 1 


State 2 


State 3 


State 1 


0.71 


0.057 


0.23 


State 2 


0.13 


0.51 


0.36 


State 3 


0.45 


0.19 


0.36 



Table 2: Mixture Means /ij/j (monthly % return) 



Component Distribution (k) 


State 1 


State 2 


State 3 


1 


1.662 


-2.588 


0.776 


2 


-30.170 


1.113 


-0.680 


Overall 


1.41 


0.097 


0.44 



Table 3: Standard Deviations ajk (monthly % return) 



Component Distribution (k) 


State 1 


State 2 


State 3 


Gaussian 1 


1.201 


2.940 


2.447 


Gaussian 2 


0.100 


2.728 


2.715 


Overall 


1.192 


2.88 


2.51 



Table 4: Initial State Probabilities (tTj) 



State 


Probability 


1 


1.0434e-024 


2 


0.99 


3 


1.6637e-017 



Table 5: Weighting Matrix (wjk) 



State 


Component k=l 


Component k=2 


1 


0.99207 


0.0079285 


2 


0.73889 


0.26111 


3 


0.7665 


0.2335 
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Table 6: Stress Test Example Scenario Tree 



Scenario Number 


Unstressed 


Stressed 


1 


-0.2521 


-10.1146 


2 


-0.0686 


-10.0433 


3 


0.0572 


-0.3384 


4 


0.8312 


1.6168 


5 


2.0156 


1.8125 


6 


2.3479 


1.8717 


7 


2.5199 


2.0075 


8 


2.6925 


2.0551 


9 


3.1681 


3.0902 


10 


3.6119 


3.0923 



Table 7: Scenario Generation Under Different Credit Risk Weightings 



Scenario Number 


Credit Risk Weighting 


0.1 


0.3 


0.5 


1 


-0.3384 


-0.3384 


-0.3384 


2 


-0.1115 


-0.1115 


-0.1115 


3 


-0.0967 


-0.1043 


-0.1043 


4 


-0.0881 


-0.0987 


-0.0987 


5 


-0.1043 


-0.0967 


-0.0983 


6 


1.6168 


-0.0881 


-0.0971 


7 


1.8125 


1.6168 


-0.0967 


8 


1.8717 


1.8717 


-0.0881 


9 


2.0075 


2.0075 


-0.0881 


10 


3.0902 


3.0902 


1.6168 
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Figure 1: FTSE-100 Empirical Monthly Returns (%) for each Month 
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From table 2's results we can see state 1 is the growth state, state 2 the recession 
state, state 3 is transitional (since the overall distribution means are positive, negative 
and about respectively). We see from table 1 that if we are in state 2, then 51% of 
the time we will remain in state 2, similarly 71% for state 1. 

We provide a GM HMM simulation of the FTSE-100 in figure 2 to compare its 
modelling performance to the empirical data (see figure 1). Simulating the GM HMM 
involves randomly selecting a Markov state in accordance to the transition matrix A 
and at each state randomly sampling a value from its pdf. Notice that the model cap- 
tures many key features of the empirical data, namely extreme losses, autoregression 
and jumps. 



Copula Approximation 

Copulas are used to see a factor's impact upon the overall multivariate distribution 
e.g. mean, variance etc... . For example, we may wish to determine the impact upon 
a FT SE-100 index tracking portfolio (or any portfolio) from credit risk or operational 
risks VoiOSj. Rather than applying copulas to analyse the impact of operational and 
credit risks, we can apply a GM model and represent these two factors as two Gaussian 
components. For example, if we choose to examine the impact on state 1 then we may 
have the distribution: 



Risk Factor 


Mean 


Sigma 


Weighting 


FTSE-100 Portfolio Return Component Gaussian 1 


1.66 


1.2 


0.69 


FTSE-100 Portfolio Return Component Gaussian 2 


-30 


0.1 


0.01 


Operational Risk 





1 


0.2 


Credit Risk 


-0.1 


0.01 


0.1 



We can also calculate the moments, the pdf is a simply weighted sum: 
fiO = 0.7iV(e, 1.41, 1.2) + 0.2iV(e, 0, 1) + 0.1iV(^, -0.1, 0.1). 



(10) 



Since the overall distribution is univariate, standard GM sampling and scenario gener- 
ation techniques can be applied. This demonstates the simplicity and computational 
advantages of GM models compared to copulas. 



Factor Specific Stress and Shock Testing (Robustness Analysis) 

Rather than conducting a worst case scenario analysis for shock testing, we can more 
accurately model a "shock event" (such as a "credit crunch" ) with a Gaussian mixture 
distribution component of chosen weighting, mean and variance. For instance we may 
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propose to stress test state I's distribution using: 

/(O = 0.9(0.99Ar(^, 1.66, 1.2) + 0.01Ar(^, -30, 0.1)) + 0.1iV(^, -10, 0.1), (11) 

where the first two Gaussians represent state I's (unstressed) distribution and the 
third represents an extreme event. Notice that the extreme event's inclusion fattens 
the left tail of the overall distribution so all scenarios sampled from this distribution 
are affected, unlike in standard worst case scenario analysis. We show a scenario tree 
generated from equation (fTTj) in table 6; notice the presence of more negative events 
scenarios due to the "shock" event Gaussian. 

Using Gaussian mixtures, we can conduct robustness analysis with respect to 1 
particular risk factor and calculate the impact on the overall distribution. For example 
we vary the credit risk's weighting by 0.1,0.3,0.5 and show the impact on scenario 
generation in table 7. As can be seen in table 7, as we increase the credit weighting 
the impact of negative scenarios becomes more prevalent. 

5. Conclusion 

In this paper, a new method of scenario generation has been proposed and demon- 
strated. We have shown how the properties of GM HMM enable us to capture impor- 
tant time varying dynamics of data, in particular reversion and jumps. This enables 
us to produce more effective scenario trees. The tractable properties of GM enable us 
to further explore SP models, through richer stress and shock testing with respect to 
particular factors. Additionally we offer a simple method of approximating copulas for 
risk analysis. We have demonstrated these properties of GM HMM scenario generation 
by numerical experiments on FTSE-100 data. 
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